## This file assesses the matches with whites using all possible DAs in canada.
options(digits = 4, width = 132)

library(here)
library(tidyverse)
library(RItools)

source(here::here("Design", "nonbimatchingfunctions.R"))

load(here::here("Design", "matches_anyDA_new.rda"), verbose = TRUE)
load(here::here("Data", "wrkdatOwnMap_anyDA_new.rda"), verbose = TRUE)

wdat0 <- wrkdatOwnMap_new %>%
  filter(!is.na(vm_change) & !is.na(social.capital01) & !is.na(community.resp01)) %>%
  droplevels()

wdat0 <- left_join(wdat0, matches_anyDA_new)
table(is.na(wdat0$pair))

## Easier to think about "higher perceiver" versus "lower perceiver" when it comes to asking whether the matching worked.
## Further simplify the working file
wdat1 <- wdat0 %>%
  filter(!is.na(pair)) %>%
  group_by(pair) %>%
  mutate(perc_more = rank(vm.community.subj) - 1) %>%
  ungroup() %>%
  mutate(pairF = factor(pair)) %>%
  select(
    perc_more, da_prop_vm_20pct_06, vm_change, pair, pairF, social.capital01, community.resp01,
    vm.community.norm2, vm.community.subj, vcid, dauid, csd_pop_06, csd_prop_vm_20pct_06,
    da_pop_dens_06, csd_urban_06, csd_pop_dens_06, community_area_km
  ) %>%
  droplevels()

xb1_ranked0 <- xBalance(perc_more ~ da_prop_vm_20pct_06 + vm_change,
  strata = list(unstr = NULL, pair = ~pair), report = "all", data = wdat1
)
xb1_ranked0$overall
xb1_ranked0$results

xb1_ranked <- balanceTest(perc_more ~ da_prop_vm_20pct_06 + vm_change + strata(pair), data = wdat1, p.adjust = "none")
xb1_ranked$overall
xb1_ranked$results[, c("Control", "Treatment", "adj.diff", "std.diff", "p"), ] # , "pair"]


## An alternative for the overall test which should be fairly close to it
library(formula.tools)
library(coin)
coin_fmla <- da_prop_vm_20pct_06 + vm_change ~  vm.community.subj | pairF
coin_test <- independence_test(coin_fmla, data = wdat1, teststat = "quadratic")
coin_test_perm <- independence_test(coin_fmla, data = wdat1, teststat = "quadratic", distribution = approximate(nresample = 1000))
coin::pvalue(coin_test)
coin::pvalue(coin_test_perm)

coin_fmla_rank <- da_prop_vm_20pct_06 + vm_change ~ perc_more | pairF
coin_test_rank <- independence_test(coin_fmla_rank, data = wdat1, teststat = "quadratic")
coin_test_rank_perm <- independence_test(coin_fmla_rank,
  data = wdat1,
  teststat = "quadratic", distribution = approximate(nresample = 1000)
)
coin::pvalue(coin_test_rank)
coin::pvalue(coin_test_rank_perm)

pairdiffs <- wdat1 %>%
  group_by(pair) %>%
  summarize(
    da_prop_vm_20pct_06 = abs(diff(da_prop_vm_20pct_06)),
    vm_change = abs(diff(vm_change)),
    vm.community.subj = abs(diff(vm.community.subj)),
    csd_pop_06 = abs(diff(csd_pop_06)),
    csd_pop_dens_06 = abs(diff(csd_pop_dens_06)),
    csd_prop_vm_20pct_06 = abs(diff(csd_prop_vm_20pct_06)),
    # da_pop_dens_06 = abs(diff(da_pop_dens_06)),
    community.area.km = abs(diff(community_area_km))
  )

sapply(pairdiffs, summary)

sapply(pairdiffs, function(x) {
  quantile(x, sort(unique(c(seq(0, 1, .1), seq(.9, 1, .01)))), na.rm = TRUE)
})

### facts for the body.tex

# number of people in final design
nrow(wdat1)
# number of DAs
length(unique(wdat1$dauid))


### The plot

## A plot of the success of the matchings to balance da_prop_vm_20pct_06 but allow variation in vm.community.norm2
pdf(file = here("Figures_Tables", "nbmplot_anyDA_new.pdf"), width = 6, height = 6)
par(oma = rep(0, 4), mgp = c(1.5, .5, 0), mar = c(3, 3, 2, 0))
nbmplot(wdat1,
  yvar = "da_prop_vm_20pct_06", xvar = "vm.community.norm2", strata = "pair", points = FALSE,
  xlab = "Perceptions of Local Community Drawings", ylab = "Census Figures",
  main = "Proportion Visible Minority"
)
abline(0, 1, col = gray(.5), lwd = 1)
dev.off()

system("touch Data/match_assess_anyDA_new.done")
